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Abstract — In this article, tlie joint fluctuations of the extreme 
eigenvalues and eigenvectors of a large dimensional sample 
covariance matrix are analyzed when the associated population 
covariance matrix is a finite-rank perturbation of the identity 
matrix, corresponding to the so-called spiked model in random 
matrix theory. The asymptotic fluctuations, as the matrix size 
grows large, are shown to be intimately linked with matrices 
from the Gaussian unitary ensemble (GUE). When the spiked 
population eigenvalues have unit multiplicity, the fluctuations 
follow a central limit theorem. This result is used to develop 
an original framework for the detection and diagnosis of local 
failures in large sensor networks, for known or unknown failure 
magnitude. 



I. Introduction 

In the field of fault detection and diagnosis, one of the 
elementary requests is the fast, reliable and computationally 
hght identification of a system failure. In dynamical scenarios, 
these systems are composed of several fluctuating parameters 
whose evolutions are tracked by a mesh of sensors reporting 
successive correlated and noisy data measurements to a central 
decision unit. With the growth in size and complexity of such 
systems, it becomes increasingly difficult for decision units 
to process simultaneously and at a low computational cost 
the augmenting load of reported measurements. Examples of 
such systems are the recent cognitive radio networks \\\ and 
smart grid technologies f2}\ . In the former, multiple cooperative 
wireless communication devices, referred to as the secondary 
network, exchange sensed data in order to decide collectively 
which communication bandwidths are left unused by the 
licensed, also called primary, network users. Fast detection 
of sudden changes, e.g. new primary user communications, is 
here demanded to minimize the interference generated by sec- 
ondary users. In the smart-grid framework, a large dimensional 
graph of interconnected electricity producers, transportation 
systems, and consumers evolve in real-time, their behaviour 
being reported by diverse sensors such as voltage phasor 
measurements |3 1 at the nodes of the electricity grid to regional 
controllers. Fast detection of link and node failures is requested 
in this scenario to minimize the risk of cascaded failures 
leading to regional blackouts |4|. There exists a rich literature 
on failure detection, diagnosis and change-point estimation, 
ranging from off-line detection methods of uncorrected data 
lEl, JH to fast change detection methods in time correlated 
signals Q, IHl, H. Subspace methods were in particular 
proposed to detect system changes from modifications in the 
eigenstructure of sampled covariance matrices for dynamical 



systems ifTol . ifTTI . In this article, we propose a novel sub- 
space approach to solve the problem of off-line detection and 
identification of local failures from independent or linearly 
time-correlated samples. 

We precisely assume the observation of measurements sug- 
gesting an error has already occurred in the network. We 
wish the detection of a failure to be fast so we will assume 
that the number n of successive sensor data reports is not 
extremely large with respect to the size N of the network. We 
will also assume that the hypothetical failure scenarios are, to 
some extent, known in advance. In this context, calling !Ho 
the hypothesis that the system does not undergo any failure 
and !Kfc, 1 < k < K, the hypothesis that a failure of type 
k occurs, the question of failure detection and localization 
consists in proceeding to the successive hypotheses tests: 
(i) decide whether the concatenation matrix of n successive 
network observations E = [si,...,s„] e £^Nxn suggests 
hypothesis CKq or its complementary J{o (i.e. the event union 
of the JCfc), and (ii) upon decision of Sq, decide what 3-Ck 
is the most likely. Both problems are optimally solved by 
multi-hypothesis Neyman-Pearson tests [12] with maximum 
likelihood performance given the observation S. However, 
this procedure is computationally intense for large system 
dimensions and large K. 

The approach under consideration here follows the theory 
of large dimensional random matrices. Precisely, we consider 
the setting where both N and n grow large and such that 
cn — N/n — > c, with < c < 1, as A^, n ^ oo. Under 
this assumption, we develop asymptotic results on the extreme 
eigenvalues and associated eigenvectors of a certain family of 
random matrices to provide novel subspace methods for failure 
detection and localization. Our interest is on random matrices 
of the spiked model type, introduced by Johnstone ifTSJI . 
specifically here of matrices modeled as S = {In + Pj^X, 
where X is a left-unitarily invariant random matrix and P 
is a rank-r Hermitian matrix with r -^ N. Such matrix 
models have been largely studied in the recent random matrix 
literature, very often in the special case where X is a standard 
Gaussian matrix, which refers in this article to a random 
matrix with independent C3\f(0, 1/n) entries. In lfT4l . ([T5 1, for 
X a standard Gaussian matrix, it is first shown that there exists 
a natural mapping between the extreme (empirical) eigenvalues 
of EE* and the (population) eigenvalues of P. It is then 
proved that, almost surely, the extreme empirical eigenvalues 
converge to deterministic limits in the asymptotic setting, 
found either at the edges of the support of the Marcenko-Pastur 



law fW\, i.e. the (almost sure) weak limit of the eigenvalue 
distribution of XX*, or away from them, depending on the 
corresponding population eigenvalues. This induces a phase 
transition having important consequences on fault detectability 
in sensor networks (TT\. This observation is extended to the 
non-Gaussian case and generalized to other spiked models in 
IITSl . The fluctuations of the extreme eigenvalues are studied 
with different approaches depending on whether the limiting 
eigenvalues are found at the edge or outside the support of 
the Marcenko-Pastur law. When at the edge, it is proved 
successively in ED, IHl, EO), lEl, |l22l that the (centered 
and scaled) limiting eigenvalue has Tracy- Widom fluctuations. 
When outside the support, those fluctuations are linked to the 
distribution of the eigenvalues of GUE matrices, as shown in 
|l21|. In the specific case where the spiked eigenvalues of P 
have unit multiplicity, the fluctuations are Gaussian. 

In this article, the properties of the extreme eigenvalues in 
a spiked model will be used to provide failure detection tests, 
in the same line as [|23l . ll24l . For failure localization, the 
information on the eigenvalue position can be used to reduce 
the number of hypotheses K. However, tests solely based on 
the limiting properties of the eigenvalues will turn out to be 
inefficient to discriminate the remaining hypotheses and we 
therefore develop novel results on the eigenspaces associated 
to these eigenvalues. In ll25l . it is shown, in the real Gaussian 
case, that the projection of the eigenvectors associated with 
the extreme empirical eigenvalues of SS* on the subspace of 
the corresponding population eigenvectors of P has a positive 
limiting norm, which is close to 1 for c small. This remark 
is extended in ifTsl to the non-Gaussian case. This property 
is the basis of our novel failure diagnosis method. However, 
the fluctuations of the eigenvector projections, fundamental 
here to derive test statistics for failure localization, have 
never been derived before for either the Gaussian or the non- 
Gaussian cases. The main mathematical result of this article. 
Theorem m provides the joint fluctuations of the eigenvalues 
and eigenspace projections for the eigenvalues found away 
from the limiting support of XX*. Our proof technique is 
largely inspired by | l8l , ll26l . lIZTl . We also use some tools 
from ll28l and ll29l . Based on these results, we suggest an 
original framework for local failure detection and identification 
in large sensor networks. 

The remainder of this document unfolds as follows. Section 
nil introduces elementary examples of sensor networks for 
which local failures translate into small rank perturbations of 
the identity matrix. Section |lll] reminds important notions of 
random matrix theory and introduces the main mathematical 
results of this article. Practical application algorithms along 
with simulations are then carried out in Section |IV] Finally, 
Section IV] concludes the article. 

Notations: In this document, capital characters stand for 
matrices while lowercase characters stand either for scalars or 
vectors, with In e i^nxn ^^^ identity matrix. The i^^ entry of 
a vector x is denoted x{i). The symbol (•)* denotes complex 
transpose. For a function / and a Hermitian matrix X E 
CA'xJV^ /(X) = C/diag(/(Ai(X)),...,/(Aw(X)))[/* with 
Ai(X), . . . , \NiX) the eigenvalues of X andU G C^x^ flie 
unitary matrix of its respective eigenvectors. The symbol St^ 



denotes the support of the probability measure tt. The notation 
Span(wi, . . . , Uk) designates the space generated by the vec- 
tors ui, . . . ,Wfc. The notation 5^ is the space orthogonal to 
S. We denote C+ = {z e C, 5(z) > 0}. The norm ||X|| of a 
Hermitian matrix X is understood as the spectral norm, and the 
norm ||a;|| of a vector x is understood as the Euclidean norm. 

as F 

The notations '—4', '=>', and ' — >' denote convergence 
almost surely, weakly, and in probability, respectively. The 
symbol E[-] denotes expectation. The notation 1^1(2;) is the 
indicator function on the set A. 



II. Detection and localization of local failures 

To motivate the study of the fluctuations of extreme eigen- 
values and eigenvectors of sample covariance matrices in the 
context of local failures in large dimensional sensor networks, 
we introduce in the following two basic examples of sensor 
network failure scenarios, which can all be modeled as small 
rank perturbations of the identity matrix, as well as related 
engineering applications. 



A. Node failure 

Consider the following model 



(1) 



where H e C^''^ is deterministic, 6 = [6{l), . . . ,6{p)Y e 
€P, w e C^ have independent and identically distributed 
(i.i.d.) complex standard Gaussian entries, and ct > 0. We 
denote y = [y{l), . . . ,y{N)]^ e C^. In a sensor network 
composed of N nodes, y represents the observation through 
the channel H of the vector 9, constituted of centered and nor- 
malized independent Gaussian system parameters!^ impaired 
by white Gaussian noise. Therefore, E[yy*] = HH*+a^lN = 
R. 

In case of failure of sensor k, y{k), the fc*'* entry of y, will 
start suddenly to return noisy outputs inconsistent with the 
model ([TJ. Assuming this noise Gaussian with zero mean and 
variance cr^ and denoting y' the observations of the network 
with failure at sensor k, we can write 

y' = {In - eke*f.)H0 + akdkeld' + (tw 

where 6' is distributed like 9 but independently and efc € C^ 
is such that efe(fc) — 1 and 6^(1) = 0, for all i ^ k. 

Therefore, y' is Gaussian (as the sum of Gaussian variables) 
with zero mean and variance 



E[2/V1 = {In - eke*k)HH*{lN ~ eke*,) + alekcl + aH^. 

Denoting s — R^^y', we have 

E[ss*] =In- R^^'HH*ekelR-^ 

+ R-hk \{elHH*ek + crl)elR-^ ~ e*,HH*R-^ 



' Up to a right-product of H by a positive diagonal matiix, the vaiiance of 
the enti'ies of 9 can be assumed all equal to one without loss of generality. 



Therefore, the population covariance matrix E[ss*] is a per- 
turbation of the identity matrix by 



Pk = R-hk 



ielHH*ek + at)elR-2-elHH*R- 



R-^HH*ekelR'' 



(2) 



Notice that the image of Pk is included in the subspace 
Spa.n{R^2ek,R^2HH*ek) and is therefore at most of di- 
mension two. Generahzing the above to M node failures at 
nodes fci, . . . , kM, the vector s is now such that 

E[ss*] ^In - R^^HH*EE*R-^ 
+ R-^E 



with E = [efcj,. 
now dU becomes 



{E*HH*E + JS?)E*R-^- - E*HH*R-2 
..,ekj,„], A = diag(crfe,,...,crfc„), where 



Pk,,...,kM - R^^-E [{E*HH*E + A^)E* - E*HH*] R-^ 
-R-iHH*EE*R-i (3) 

which is a matrix of rank at most 2M. 

B. Sudden parameter change 

Consider again the elementary model ([TJ and now assume 
that, instead of a sensor failing, d{k), the k*^ entry of 6, 
experiences a sudden change in mean and variance. The 
resulting observation y' can be modeled as 

y' = H{Ip + akCkcDe + ^ikHck + (TW 

for some real parameters iik^ ctk, and where ek G C^ is defined 
by efe(fc) = 1 and ek{i) = 0, i ^ k. For this particular model, 
we may or may not suppose that Hk and ak are a priori known 
to the experimenter. In this scenario, we now have that y' is 
Gaussian with zero mean and variance 

E[yV1 = H{Ip + [ill + (1 + akf - lWel)H* + a^I^. 

Denoting R = HH* + (j'^Im as in the previous example and 
taking s = R^2y\ we finally have 

E[ss*] =In + b4 + (1 + "fe)^ - l]R^^HekelH*R-^2 

which is a rank-1 perturbation of the identity matrix by the 
matrix 

Pk^PkR-^HekelH*R-i 

with Pk ~ nl + (1 + Qffc)^ — 1. Note that, in this scenario, 
the eigenvector of Pk associated with the non-zero eigenvalue 
is independent of fik and ak- For practical applications, this 
has the interesting advantage that simple localization can be 
performed even if fik and ak are unknown. This is further 
discussed in Section ITVl 

The derivation above generalizes to sudden changes of 
multiple parameters. If the means and variances for the sen- 
sors ki, . . . , kM are modified simultaneously with respective 
parameters ^fe^ , . . . , fikM and ak,,---, ak^^ then 



E[ss*] 



LN 



R-'!HEAE*H*R-'! 



with E = [ck, ,..., ekj,,] and A = diag{(3k, ,..., PkM )- hi = 
/^fc + (1 + "^fci)^ ~ 1> which is a rank-Af perturbation of the 
identity matrix by the matrix 

Pku....kM =R'^HEAE*H*R-i. 

Note that, contrary to the one-dimensional case, the eigenvec- 
tors of Pki,...,kM depend here explicitly on the parameters ^ki 
and aki- 

In the following section, we introduce the novel detection 
and localization framework and we discuss engineering appli- 
cations referencing the examples described in this section. 

C. Detection and localization in sensor networks 

For either of the models above, let us assume a gen- 
eral scenario with K possible failure events, indexed by 
1 < k < K and let now si,...,s„ be n successive 
independent observations of the random variable s. We denote 
S = -^[si, . . . , s„] G C^^". From the fact that s is Gaussian 
with zero mean and covariance [I^ + Pk) for a certain fc, we 
can write 

^^{lN+Pk)--X 

where X E (j^wxn j^ ^ given matrix with independent 
Gaussian entries of zero mean and variance 1/n. We also 
denote for simplicity Po = for the extra event fc = 
corresponding to the no-failure scenario. 

The natural approach to detect and identify a failure event 
in a sensor network upon the observations si, . . . , s„ is to 
systematically perform a maximum likelihood test on the 
K +1 hypotheses ^q, . . . , 3-Ck, with !Kfe defined as the event 
s ~ C3\f(0, Ij\[+Pk)- However, this optimal approach has some 
intrinsic limitations. From a computational aspect, evaluating 
the probability of each hypothesis k requires to evaluate the 
term tr 'S*{In + AO ^S, an operation whose cost is of order 
0{N^) (which can be brought down to 0(A^^) using matrix 
inversion lemmas). When the number of hypotheses K and the 
system size N are large, these operations become extremely 
demanding. Pre-calculus of the inverses {In + Pk)~^ also 
requires possibly large memory storage. 

Since the node failure information is entirely captured by 
the perturbation matrix Pk, we provide in the following a 
suboptimal test relying on the properties linking Pk to the 
observation matrix S, for large system dimensions {N,n). 
Precisely, based on recent advances in the field of large 
dimensional random matrix theory ll30l . we provide a two- 
step approach to successively (i) decide on the existence of a 
failure from the location of the extreme eigenvalues of SS* 
and (ii) identify the failure event from eigenspace projections. 
This diagnosis framework relies on the asymptotic statistics 
of these extreme eigenvalues and eigenspace projections. This 
subspace approach has multiple advantages compared to the 
optimal hypothesis testing method discussed above. From 
a computational aspect, step (i) requires to determine the 
eigenvalues of EE*, hence a singular value decomposition. 
This step already provides sufficient information for step (ii) to 
become computationally cheap: on the one hand, the position 
of the extreme eigenvalues of SS* may be used to reduce 
the set Oil , . . . , CKk to a possibly small subset of consistent 



hypotheses; on the other hand, for the remaining hypotheses, 
the localization test will merely consists in the characterization 
of eigenvector projections, an operation of computational cost 
0{N). No matrix inverse needs to be computed and only the 
eigenvectors and non-zero eigenvalues of Pk need to be stored. 
The technique also has the advantage to be consistent in its 
usage of eigenvalues and eigenspace projections to perform 
hypothesis tests. Finally, as will be discussed in Section |IV] 
the framework can be extended to account easily for unknown 
failure amplitudes, which would be much more involved from 
a maximum-likelihood approach. 

Before presenting our main results, we briefly discuss 
practical applications in the failure diagnosis context. The 
application of the aforementioned results to failure diagnosis 
in sensor networks consists in deriving subspace methods for 
failure detection and identification when the system undergoes 
sudden changes that can be modeled as small rank perturba- 
tions of the data covariance matrix. The scenario of Section 
III-AI mav be exploited in particular to rapidly detect failures 
in sensors that may suddenly return inconsistent data. Fast 
detections of such sensor errors are of fundamental importance 
when decisions on system actuators are taken from these data. 
The scenario of Section ITl-B I finds applications in the cognitive 
radio setting where Gaussian baseband signals 9{1), . . . ,9{p) 
arising from a set of p fixed access points, forming the so- 
called primary network, are sensed by a network composed of 
N cooperating sensors, constituting the secondary network. 
Assuming that at a reference time instant a subset T of the 
p access points transmit data, the communication channel 
H e C^^^ between primary and secondary networks (which 
can be assumed static during a channel coherence interval) 
will contain non-zero entries in the columns indexed by CP. 
The objective of the secondary network is here to rapidly 
detect and identify changes in H (column elimination or 
creation), corresponding to evolutions of the subset CP of 
transmitting access points. These fast change diagnoses will be 
used to exploit the available radio spectrum as well as to avoid 
secondary transmissions to interfere primary communications. 
We also mention that in [31 1, a similar system model for failure 
diagnosis in power systems is analyzed. In this setting, the 
observation vector y is composed of voltage measurements 
at the N network nodes, the variable 9 contains the current 
inputs at the nodes which naturally vary due to the presence of 
unreliable renewable energy production units, and the transfer 
matrix H is closely related to the square matrix containing at 
entry {i,j) the inverse impedance of the power line connecting 
i to j (and equal to zero if i and j are not connected). The 
objective in this system is to detect local power line outages 
corresponding to impedance changes in H, which can be seen 
as rank-two perturbations of the population covariance matrix. 

The following section is dedicated to the study of the 
asymptotic eigenvalue and eigenspace projection statistics as 
the dimensions of the matrix S grow large. 

III. Main results 



borrow some of the arguments of 
to ours. 



whose context is close 



A. Notations, assumptions and basic results 

We start by summarizing the major notations and facts 
needed here. We consider a generic small rank perturbation 
model and define 

with X e C^^" left-unitarily invariant, and where the rank-r 
Hermitian matrix P has the spectral factorization P = UilU* 
with 



n 



Ujjjt 



JtA 



> Ws > > Ws+i > . . . > Wf > —1. Of course, 

: r. We write accordingly U = \Ui ■ ■ ■ Ut\ where 

We denote by Ai > • • • > AJ^ 



and oji > . . 

JlH hjt = 

Ui E C^^^' . We denote by Ai > • • • > Ajv the eigenvalues 

of SS*. For i e {1, . . . , s}, we let %{i) = ji H h j^-l, 

taking by convention jo = 0. For i £ {s + 1, . . . ,t + 1}, we 
let 3C(i) — N — {ji + ■ ■ ■ + it). One of the purposes of this 
section is to establish an asymptotic relation between LOi and 
the \-x{i)+e. for t = 1, . . . ,ji which holds under a condition on 
uJi that will be specified. We also denote by 11^ the orthogonal 
projection matrix, when it exists, on the eigenspace of SS* 
associated with the eigenvalues {^x(i)+£}'i^i- Similarly, we 
denote by Hi = UiU* the orthogonal projection matrix on the 
eigenspace of P associated with the eigenvalue uji. Finally, 
we denote by Q{z) = {XX* — zl^)^^ the resolvent of the 
matrix XX* and by a{z) = j^ii Q{z) its normalized trace, 
both analytical on C+. 

In the remainder of the paper, we shall consider the asymp- 
totic regime where n — )■ cx) and N/n — )• c € (0, 1). The 
notation n — > oo will henceforth refer to this asymptotic 
regime. 

We now state our basic assumptions: 

Al The probability law of X is invariant by left multi- 
plication by a deterministic unitary matrix. 
Thanks to the left unitary invariance of X, Q{z) writes 
as Q{z) — W{K — zIn)^^W* where A is the matrix of 
eigenvalues of XX*, W is a unitary random matrix Haar 
distributed on its unitary group, and W and A are independent. 

A2 For every z e C+, a{z) converges almost surely to 
a deterministic function m(z) which is the Stieltjes 
transfornu of a probability measure n with support 
[a,b] C (0,cx3). 

A3 We have II a: a:* 1 1 ^ 6 and (||(A:a:*)-1||)-i ^ a. 
This last assumption implies in particular that A2 is satisfied 
for all ze C\ [a,b]. 

The most classical model of a matrix X that satisfies Al- 
A3 is when X is standard Gaussian, i.e. with independent 
CA^(0, 1/n) elements, as introduced in the system models of 
Section III-AI and Section III-BI For this model, the limiting 



The derivation arguments found in this section follow the 2^6 recall that the Stieltjes transform miz) of a real measure tt is defined 



ideas of 



, and ll27l . In our proofs, we shall also for z outside the support of vr by m,(2) = / 



:dM\). 



probability distribution n is the well known Marcenko-Pastur 
distribution ||T6l . Its Stieltjes transform m(z) is given by 



m{z) 



2zc \ 



1 



V^ 



z)2-4zc (4) 



where the branch of the square root is chosen such that 
77i(C"'") C C+ and m is analytic on C \ [a,b], where 
a = (1 - Vc)^ and 6 = (1 + ^/cf. 

The unitary invariance of X is the basis of the following 
important lemma, shown in \29) using an inequality of BSI 
which involves Haar unitary matrices: 

Lemma 1: Assume Al. Let u E C^ and v £ C^ be two 
vectors with norm ||ii|| = ||u|| = 1. Denote by a{XX*) the 
eigenvalue spectrum of XX*. Given e > and z e C \ [a — 
e,b + e], denote by dz the distance from z to [a — e,b + e]. 
Then for any p > 0, 



expect the solutions of the equation det_ff(a;) — which 
are outside [a, b] to coincide with the limits of the isolated 
eigenvalues of SS*. 

Let us now study the behavior of the solutions of this 
equation. Let h{x) — xm{x) on M \ [a,b]. Since a > 0, we 
have 

The function h{x) is therefore increasing on M\ [a, b] and with 
limit as a; — )■ and — 1 as x ^ oo. Therefore, for w^ > 0, 
(01 leads to 

h{p) + ^±^=0 (6) 



h' (x) = {xm{x)) 



E [\la{XX*)cla-e,b+e]U* {Q{z) - a(z)/jv)w| 



< 



K„ 



dlNP/^ 



where the constant Kp depends on p only. Similarly, for any 

z, z' G C \ [a — £, 6 + e], we have 



E 



< 



lcr(XX*)c[a-e,b+e]"* 



o,.)0(=',.i:0(|M,„ 



p-\ 



having a unique real solution pi satisfying pi > b if and only 
if /i(5+) + (1 + uJi)/uJi < 00 When w^ < 0, Q has a unique 
solution < Pi < a if and only if h{a^) + (1 + oji)/uji > 0. 
We therefore have the following result, for which a rigorous 
proof is found in f32l : 

Theorem 1: Assume A1-A3. Let p be zero or the maximum 
index such that cjp > and h{b^) 



(1 + U}p)/Ujp 



< 0. For 



d^d^NP/^' 



i = I, . . . ,p, let Pi be the unique solution of (|6]l such that 
Pi > b. Then, 

for i = 1, . . . ,p and £ = 1, . 



■,Jt 



while 



We now start our analysis of the extreme eigenvalues and 
eigenspace projections of SS* by studying the first order 
behavior 

B. First order behavior 

1) Eigenvalues: Suppose that a; G M is not an eigenvalue 
of XX*. We first write 

det(SS* - xIn) 

= det(/jv + P) det(XX* - xIn + x[In - {In + P)~^]) 

= det{lN + P) det{XX* - xIn) 

X det(/Ar + xP{In + Py^XX* - xInT^) 

after noticing that In - {In + PY^ = P{In + PY^- 
Therefore, if x is an eigenvalue of SS* but not of XX* , it 
must cancel the rightmost determinant. This determinant can 
be further rewritten 

det(/Ar + xP{In + Py^Q{x)) 

= Aet{Ir + xVLU*{lN + UVLU*)-'^Q{x)U). 

From the identity U*{In + UVtU*)-^ = {U + VL)-^U* , we 
then have 

det(EE* - xIn) = det(/jv +P) Aei{XX* -xIn) det(^(x)) 

where H{z) = /,. + zO(/^ + Vl)-^U*Q{z)U . 
Given any z G C \ [a, 6], Assumptions A1-A3 in conjunction 
with Lemma[T]show that H{z) is defined for all n large, almost 
surely, and converges with probability one to 



A 



3<:(p+i)+i 



Let g be t + 1 or the minimum index such that Wg < and 
h{a^) + {l+ujq)/ujq > 0. For i — q,...,t, let pi be the 
unique solution of ^ such that pi < a. Then, 



for i = q, . . . ,t and £ = 1, . . . ,ji 



while 



A 



■K(q) 



In the remainder of the article, the variables wi , . . . , Wp and 
ujq, . . . ,ujt satisfying the conditions of Theorem [1] will be said 
to satisfy the separation condition. 

When X is standard Gaussian, applying Theorem [T] shows 
after some simple derivations the following result: 

Corollary 1: Consider the setting of Theorem [T] Assume 
additionally that X is standard Gaussian. Let p be zero or the 
maximum index for which ujp > y/c and q he t + I or the 
minimum index such that ojq < —y/c. Then 

1 + UJi 



A 



for i G {1, . 



DC{i)+e 

p,q,-- 



-^^ Pi = i- 

.,t),£=l,. 



(7) 



A 



3C(p+l) + l 

A, 



H{z) = Ir + Zm{z)n{Ir + n)-^ 



(5) 



(take u and v in Lemma [T] as any couple of columns of U , 
take p > 2 and use Borel-Cantelli's lemma ll33l ). We therefore 



bJi 
.,ji, while 

{i+v^r 

\3C(q) ^ (1 - Vc)^ 

Corollary [T] implies that, for Ui sufficiently far from zero 
(either positive or negative) or, equivalently, for c sufficiently 
small, the spectrum of EE* exhibits ji eigenvalues outside the 
support S*^ of the Marcenko-Pastur law tt which all converge 

^We denote by 3S'+ and x^ any quantity infinitesimally greater and smaller 
than the real x, respectively. 



to Pi. For failure detection purposes, upon observation of 
E, we may then test the null hypothesis Y, ~ X (call it 
hypothesis JCq) against the hypothesis E — (/jv + P)^ X 
(call it hypothesis Jio), depending on whether eigenvalues of 
EE* are found outside St^. Depending on the scenario, for c 
small enough, it may be that a mere evaluation of the number 
of eigenvalues outside the support suggests the number of 
simultaneous failures in the sensor network. This is the case of 
the two failure scenarios described in Section ITl- Al and Section 
III-BI However, the information on the extreme eigenvalues of 
EE*, if sufficient for failure detection purposes, is usually not 
good enough to perform accurate failure localization. This is 
because different failure scenarios, characterized by different 
perturbation matrices P, may exhibit very similar eigenvalues. 
Also, if the failure amplitude is a priori unknown, then 
eigenvalues are in general irrelevant to discriminate between 
failure hypotheses; see the application Section IIV-CI In such 
scenarios, we then need to consider eigenspace properties of 
P. This is the target of the following section. 

2) Projections on eigenspaces: Given i < t, we now 
assume that Ui satisfies the separation condition. Given two 
deterministic vectors foi , 62 G C^ with bounded Euclidean 
norms, our purpose is to study the asymptotic behavior of 
65; Hi 62. We shall show that this bilinear form is simply related 
with blllib2 in the asymptotic regime. 

Our starting point is to express blllib2 as a Cauchy integral 
||34I . Denoting Gi a positively oriented contour encompassing 
only the eigenvalues Xjc(i)+i of EE* for £ = 1,. . . ,ji, we 
have 

6^11,62 



= (b bUT^T.* - zlNy^b2 dz 

X [XX* ^ zIn + zP{In + Py^r^lN + Pyh2 dz. 
Using Woodbury's matrix identity, we have 

[XX* - zIn + zP{In + P)-Y^ 

= Q{z) - zQ{z)U [Ir + zn{Ir + n^^U* Q{z)U] "^ 

X n{ir + n)-^u*Q{z) 

= Q{z) - zQ{z)UH{z)-^n{Ir + n)-^U*Q{z) 
and taking 

ai{z)* ^zbl{lN+Pr^Q{z)U 

a2{z) = ii{ir + n)-^u*Q{z){iN + pyh2 

we obtain 

&in,&2 

= -^^j> bl{lN + P)-^Q{z){lN + P)-h2dz 
1 



2Tn 



di{z)*H{z) ^d2{z) dz. 



(8) 



By Assumption A3 and Theorem [T] with probability one for 
all large n, the first term on the right-hand side is zero, while 
the second is equal to 

(B ai(z)*H(z)^^a2(z)dz 

where 7j is a deterministic positively oriented circle away from 
[a, b] enclosing pi but none of the pj, j ^ i. Using Lemma [T] 
in conjunction with the analyticity properties of the integrand, 
one can show that ai{z)* H{z)~^d2{z) converges uniformly 
to ai{z)* H{z)^^a2{z) on 7^ in the almost sure sense, where 

ai(z)* = zm{z)bl{lN + Py^U 
a2{z) = m{z)Vt{Ir + n)-^U*{lN + P)"H2. 



0, where 

-1 



a2(z)dz. 



It results that 5^11^52 - Ti 

T^^^i ai{z)*H{z) 

Details can be found in ||29I in a similar situation. Let us find 
the expression of Ti. Noticing that 



H{zr 



t 
£=1 



1 



zm{z)-^_ 



(9) 



where 



Jji+.-.+jf-i 



Jjf+1+...+jt 



eC 



we obtain 

T, 



t 



UJl 



tt(l+^^) 



b*Ufh 



1 



2"1 



£02 



z'm?(z) 






blUeb2 1 



27ri .1^, 1 
zm'^{z) 



^. -^ , zm(z)-ri 



-dz 






zm{z) 



-dz. 



From the discussion prior to Theorem [T] it is clear that the 
denominator of the integrand has a zero in the interior Int(7i) 
of the disk delineated by ji only when i = i, and then this 
zero is simple due to the fact that h'{x) ~ {xm{x))' never 
vanishes on R \ [a, b]. Applying the residue theorem [i34j and 
observing from ^ that 1/(1 + w^) = (1 + h{pi))/h{pi), we 
obtain the following limits: 

Theorem 2: Assume A1-A3. Given i < t, assume that tUi 
satisfies the separation condition. Let &i e C^ and 62 G C^ 
be two sequences of increasing size deterministic vectors with 
bounded Euclidean norms. Then 



where 



bln,b2 - C^K^^b2 ^ 

_ m{p,){l + h{p,)) 







h'iPi) 

In particular, we find after some derivations: 
Corollary 2: Under the assumptions of Theorem |2l let X 
be standard Gaussian. Then 

.-2 



blU,b2 



1 



1 



cw, 



^blU,b2 ^^ 0. 



This result is consistent with ||25l derived in the real Gaussian 
case for eigenvalues with unit multiplicity. 

Theorem |4] and Corollary |4] provide an interesting charac- 
terization of the eigenspaces of P through limiting projections 
in the large dimensional setting. In the context of local failure 
in large sensor networks, it is therefore possible to detect and 
diagnose one or multiple failures by comparing eigenspace 
projection patterns associated with each failure type. Precisely, 
an appropriate diagnosis consists in determining the most 
likely failure type among all hypothetical failures, given the 
extreme eigenvalues and associated eigenspace projections of 
SS*. To this end though, not only first order limits but also 
second order behaviour need be characterized precisely. This 
is the target of the following section. 

C. Second order behavior 

Before studying the fluctuations of \jc(i)+e, i = 1, ■ • ■ , ji, 
when uji satisfies the separation condition, we first remind 
for later use in our applicative framework the fluctuations of 
^x(i)+i when uJi does not satisfy the separation condition, and 
when X is a standard Gaussian matrix. For this, we have the 
following theorem [20.1, [22.1, [21]. 

Theorem 3: Let X be standard Gaussian, then if < w, < 



7V3. 



A 



3C(i)+i ■ 



(l + Vc)' 



{l + Vc)^y/c 



and, if — \/c < cji < 0, 



iV3 



>^x{i}+i - (1 - Vcf 



To 



To 



-(i-VS)M 

for i! = 1, . . . , ji, as n — > cx), where T2 is the complex Tracy- 
Widom distribution function |fT9l . 

The tools used to derive Theorem[3]are much different from 
those exploited here and will not be discussed. Note that in 
||35l , an extension to the case where X may be correlated 
is provided but only considers the fluctuations of the largest 
eigenvalue. Similar to [24|, Theorem [3] will be used to derive 
tests to decide on the presence of eigenvalues outside the 
support of the Marcenko-Pastur law. For failure detection 
purposes in sensor networks, this will be used to declare a 
failure prior to diagnose the fault. Then, to diagnose a failure, 
second order statistics of both eigenvalue and eigenspace 
projections when the separation property arises are needed. 
This is the aim of the remainder of the section. 

We now turn to the second order analysis of the eigenspec- 
trum of SS* when LOi satisfies the separation condition, and 
when X is only assumed to satisfy A1-A3. We first need the 
following additional assumption: 

A4 For an z e C \ [a, 6], 

TD> 





N {a{z) ~ m{z)) 



In practice, this assumption means that the fluctuations of 
the spectral measure of XX* are negligible with respect 
to those of the A^ and h\Iiib2, which will be shown to 



be of order V^- This assumption is satisfied by most of 
the random matrix models of practical importance in our 
context, provided \/N{N/n — c) — > 0. The classical illus- 
trating example in this regard concerns the standard Gaussian 
case. Denote by m„(z) the Stieltjes transform (|4|i where c 
is replaced with N /n, and let 7r„ the Marcenko-Pastur law 
associated with r7i„(z). For any z e C \ [a, &], the function 
f{x) = [x — z)"^ is analytic outside the support of 7r„ for 
n sufficiently large. As a consequence. Theorem 1.1 of ll36ll 
shows that \/ N{a{z) — m„(z)) — > 0. Assuming in addition 
that y/N{N/n — c) -> 0, it is not difficult to show from ^ 
that ^/N{mr^{z) - m{z)) -^ 00 

For practical purposes, we shall also assume: 

A5 Each w^, 1 < i < i, satisfies the separation condition. 

The main result of this section is the following theorem. 

Theorem 4: Assume A1-A5. For i = 1, . . . , i, let 



V^, 



Nu* (n, - qin) u. 



and 



A 



3C(i) + l 



^{i)+ji 



For p e M \ [a, b], let 
Dip) 



h{p)ii+h{p))h" jp) 

P 

h'ip) 



hip){i+h{p)y 




and 



Rip) 



m'(p) — m{p)^ 
m"{p)/2 — m(p)m' (p) 



m" {p)/2 — m{p)m' [p) 
m^^'>{p)/Q-m'{pf 



where m^'^' is the third derivative of m. Consider the matrices 



({D{p,)R{pi)D{p,rY®^ 



Ml., 
M2,, 



where A/1,1, M2.1, ■ ■ ■ , Mi^t, M2,t are independent GUE ma- 
trices such that Mii and M2.i are both ji x ji matricesO Let 
Li be the ffi.-'' -valued vector of eigenvalues of Ki arranged in 
decreasing order. TherO 



\(yi,n, Li,n))i 



{{G^,U))\ 



(10) 



Remark 1: To be more precise, (Vi_„, i^ „) can be cast into 
an Wi ^^^ -valued random vector after rearranging the real and 
imaginary parts of the elements of Vi^n and taking into account 
the Hermitian symmetry constraint. Hence, the convergence in 
law specified by (fTOl i acts on the space of probability measures 
on R'^+Ji+-+J*. 

It is useful to note that when the convergence rate of N/n is slower than 
1/yN, our analysis remains true if we replace m{z) with a finite horizon 
deterministic equivalent 1 37 1, 1 38 1. For simplicity, we have chosen not to enter 
into these details here. 

^We remind that matrices from the Gaussian unitary ensemble, or GUE 
matrices, are random Hermitian matrices with independent standard real 
Gaussian diagonal entries and independent standard complex Gaussian upper- 
diagonal entries |39|. 

*It is clear that if AS is not fulfilled. Theorem |4] generalizes by letting the 
index i in UOt span only the subset {1, . . . , p, g, . . . , t}. 
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Fig. 1. Empirical and theoretical distribution of the fluctuations of mi with 
r = 1, X has i.i.d. zero mean variance 1/n entries, N/n = 1/8, N = 256 
and cJi = 1 (top) or lji = 0.5 (bottom). 



Theorem |4] provides a very general expression of the joint 
limiting fluctuations of both eigenvalues and eigenspace pro- 
jections. It is particularly interesting to note that the fluctua- 
tions of {Vi,n,Li^n) are asymptotically independent across i. 

This theorem is an interesting new result to the field of 
large dimensional random matrix theory, although it might be 
difficult to use in practice since one needs to derive explicitly 
the joint density of {Gi,Li)l^^. Nonetheless, Theorem |4] can 
be immediately put to practice in two scenarios. The first 
scenario corresponds to the case where only the fluctuations 
of the eigenspace projection vector (l^,n)i=i is of interest. In 
this case, {Gi)l^i is a correlated Gaussian random variable. 
Reminding that Mi,i can be seen as an R-'i' -valued Gaussian 
vector with entries of zero mean and unit variance, (Gi)l^^ 
can then be seen as a random Gaussian vector in RJi +■■■+.?* 
with entries {jf_i — 1) to jf of zero mean and variance 
[D{p,)R{p,)D{pi)*]j^^, for each i e {1, . . . ,0 wifli jo = 0. 
The second scenario of interest, which is discussed at length in 
the following, corresponds to the case where the multiplicities 
of the eigenvalues of P are all equal to one. In this case. 



Mii and M2.i are independent Gaussian variables and we 
immediately have the following corollary: 

Corollary 3: Consider the setting of Theorem|4] Assume in 
addition that ji = 1 for all i. Then 



with 



R 



{{Vi,n, -^i,n))j^i 



D{pi)R{p,)D{piy 



?^(0, R) 



D{pr)R{pr)D{pr)* 

After some calculus, in the standard Gaussian case, we 
further have: 

Corollary 4: Under the assumptions of Corollary [3] if ^ is 
a standard Gaussian matrix, then D{pi)R{pi)D{pi)* ~ C{pi), 
where 



C(P.) ^ 



' (1+' 






1 



iTcl 



{uJi+c)^ZJi 



(l+a^,)^c^ 



Due to its simple expression. Corollary |4] is particularly 
handy to use in the context of failure diagnosis when hypo- 
thetical failures are characterized by distinct values of LOi, as 
will be shown in Section HV] 

In Figure [T] the histogram of a simulation of 10 000 re- 
alizations of the projection V^i_„ — ^fN{\u\ui\'^ — (,i), with 
ui — Ui ^ C^, uiu{ = Hi of unit rank, and X standard 
Gaussian, is depicted against the asymptotic Gaussian law 
derived in Coroflary H for r = 1, N/n = 1/8, N = 256 
and oji successively equal to 1 and 0.5. In this scenario, 
\J N jn ~ 0.35. For wi = 1, the simulation shows a rather 
accurate fit between asymptotic theory and simulation. For 
wi — 0.5, the Gaussian approximation is much less accurate. 
This is due to the value uj\ = 0.5, which is rather close to 
^J N/n ~ 0.35. The value N = 256 is here insufficient for the 
large dimensional behaviour of the fluctuations of |u5;uip to 
arise. This behaviour will have important consequences for the 
question of diagnosing failures which are difficult to observe. 

The remainder of this section is devoted to the proof of 
Theorem |4] We start with the following lemma, which deals 
with the asymptotic behavior of the Vi^„. This lemma will be 
proved in Appendix lA-AI 

Lemma 2: Assume A1-A5. Let 



V,n = 



h{p,){i + h{p,)) 



m{pi)I)Ui 



Then for any i e {1, 



U*{Q'{p/)~m'{pi)I)U, 



..,t}. 



V,.. 



V^, 



0. 



We now consider the isolated extreme eigenvalues. In order 
to study the asymptotic behavior of these eigenvalues, we shall 
adapt to our situation the approach of ||26| . For i = 1 , . . . , i, 
consider real numbers xi{i) > yi{i) > X2{i) > J/2(i) > 
■> Ujii'i')- Since the separation condition is satisfied by 
assumption for each i — 1, . . . ,t, the equation dct H{x) ~ 



has r roots outside [a, b] with probabihty one for all n large. 
Therefore, we have the equivalence relation 



^xi{i) > VN{Xjc[.t)+i - Pi) > ytH), 

I — i, . . . , t, t — i, . . . , J, 



^ 



Nn aetH[p, + ^) det h(p, + ^ ) < 0, 



N 
i^l,...,t, i^l 



N 

I ■ • • I Ji 



(11) 



This equivalence leads us to study the limits of the finite 
dimensional distributions of the R* -valued random process 

N^'^^ detH{pi + x/VN) in the parameter x. This is 

given by the following lemma, proved in Appendix lA-BI 

Lemma 3: Assume A1-A5. Define the ]R*-valued random 
process 

Xn{x) = 

£^i 

X det (^VNp,p,U*iQ{p,) - m{p,)I)U, + /3,/i'(p,)x/)] _ 
where f3i — uJi/{l + lui). Then 

{Xn{xi),...,Xn{Xp)) > 

■ : Xn). 



V 



NU*{Q{p^)-m{p,)I)U, 



for every finite sequence (a:i, . 

Let -B be a rectangle of R^i^ ^^t . The real and imaginary 
parts of the elements of the Vi „ defined in Theorem |4] can 
be stacked into a M-'i^ '"^* -valued vector Vn when taking 
into account the Hermitian symmetry constraint. A vector Vn 
with the same size can be constructed similarly from Vi_„ 
defined in Lemma |2] Let i„ be the R''-valued vector L„ = 
[L]^„, . . . , LJ^ (see Theorem lU and let C be the rectangle 
of R*" determined by the left hand side of ( fTTT i. so that this 
event is written [L„ G C]. Let Z„ be the vector obtained by 
arranging the eigenvalues of the matrices 

T - P» 

for i = 1, . . . , t, similarly to the elements of i„. From Lemma 
12] Lemma |3] and the discussion preceding Lemma [3] we have 

P[14 e B,i„ e C] - P [Vn e B,Z„ e C] ^ 

as n — > oo, for arbitrary rectangles B and arbitrary rectangles 
C specified at the left hand side of ( fTTT i. Observe that 

Ti,n\ ■' \U*{Q'{pi)-m'{pi)I)Ui 

In order to terminate the proof of Theorem|4] we shall make 
use of the following lemma, that we state in a slightly more 
general form than needed here. 

Lemma 4: Assume A1-A4. Let fi, ■ ■ ■ , ft and gi, . . . ,gthe 
real functions analytical on a neighborhood of [a,b]. Let 5„ 
be the t-uple of random matrices 



Sn = Vn 



u*Mxx*)u,-(jMX)d7T{\))L 



U:g,{XX*)U,-{Jg,{X)d7r{X))L 



For i 



1, . . . ,t, define the covariance matrices 

[{ft-JftdT^) {9i-j9idTr)] 



gi - J gidn 

Then Sn converges in distribution towards 

t 



dn. 



{Rf(^ij,) 



M2, 



(12) 

1=1 

where the matrices (Afi.i, Af2,i, • ■ • , Mi,t, Af2,t) are indepen- 
dent GUE matrices such that A/i ^ and A/2,i have dimensions 

ji X ji- 

The proof is provided in Appendix lA-CI 

Applying this lemma with fi{X) = 1/(A — pi) and gi{X) = 
1 / {X— pi)"^ , Ri takes the value R{pi) provided in the statement 
of Theorem |4] It results that 



Vi_ 



i=l 



i=l 



where the convergence takes place on the space of probability 
measures on M^'^i^ ^^t\ This completes the proof of The- 
orem |4] 

IV. Application 

In this section, we provide a general framework for local 
failure detection and diagnosis in large sensor networks, such 
as the examples proposed in Section |II] based on the results 
of Section Hn] This framework is a two-step approach for 
successively (i) detecting failures within a given maximally 
acceptable false alarm rate and (ii) upon positive detection, 
diagnosing the failures with high probability. Simulations are 
then run to validate the proposed algorithms. 

We recall that we assume a number K of failure scenarios 
indexed by 1 < fc < A'. Scenario k is characterized by the 
population covariance matrix In +Pk. We additionally denote 
Pq — for convenience. We consider precisely the detection 
and localization to be made for failure models of the type 

S = ilN+Pk)^X, l<k<K, with Pk = E-il OJk,rUk,rUl. 

of rank r*. = Y.lLiJk,i, where UkA e C^""^"* and ujkA > 
... > Uk,sk > > iOk,sk+i > ... > uJk,tk, and X standard 
Gaussian. The eigenvalues of EE* are denoted and ordered as 
before by Ai > . . . > Xn- We denote J-Cq the null hypothesis 
for the model 'S — X and JC^ the hypothesis for a failure 
of type k. For the failure scenario k, we also take pk to be 
zero or the smallest index i such that ujk,i > \/c, and qk to 
be tk + 1 or the largest index i such that Wk.i < —\fc. 

Remark 2: This model (and in particular the models of 
Section lTl-Al and Section lTl-Bb may be extended by introducing 
deterministic linear time correlation between the successive 
observations si,...,s„. That is, we can write X = XT 2 
for some time-correlation matrix T and a standard Gaussian 
matrix X. In this scenario, X being left-unitarily invariant, 
the localization scheme proposed extends naturally to this 
scenario, as long as the eigenvalue distribution of T converges 
weakly to a compactly supported distribution as N grows 
large and that T has asymptotically no eigenvalue outside this 
support. Alternatively, a whitening procedure may be used 
prior to failure detection using S' = ST^^ as the random 
matrix under study, provided that T is invertible. 
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A. Detection algorithm 

As stated in the introduction, the detection phase relies on 
existing resuhs, and more specifically on the fluctuations of 
the largest eigenvalues given by Theorem [3] The detection 
algorithms proposed here parallel that introduced in [24| in 
the context of collaborative signal sensing. The objective is to 
decide between hypothesis JCq and its complementary JCq. 

First assume that all Pk only have non-negative eigenvalues. 
From Theorem [1] the largest eigenvalue Ai of SS* tends to 
the right edge of the support S"^ of the Marcenko-Pastur law 
for all large n under S^q, while Ai is found away from this 
edge under 'Kq ;/ the largest eigenvalue ujk.i of Pk exceeds 
\/c. We will therefore assume in the following that Wfc i > \/c 
is verified for all k. That is, we assume that ci^ ~ N/n < c+ 
where c+ is defined as 



c+ 



= mi{ujli, l<k<K}. 



This condition allows for a theoretically almost sure error 
detection, as N,n ^ oo. We then rely on Theorem [3] to design 
an appropriate hypothesis test. Our test consists in rejecting 
hypothesis J-Cq if the probability in favor of JCq is sufficiently 
low. That is, for a given acceptable false alarm rate ryO the 
statistical test is defined as 



IKq 



J-Cq 



(13) 



where X[ is given by 



A'l ^ TVi 



Ai-(1 + V5^)^ 

{1 + y^)^C% 



That is, the test verifies whether Ai exceeds some threshold 
above which the probability for 5{q is less than t]. 

If the matrices P^ are now all non-positive definite, then, 
symmetrically, we need to set cjy such that the smallest 
eigenvalue A at of SS* is visible on the left-hand side of St^. 
That is, we take N, n to be such that cat < c_, with c_ defined 
as 

c_ =inf {wLfc> 1 < fc </-i:}. 



The decision test is in that case given by 



A^ ^ {T2)-'{l-v) 



J-Cq 



where A^ is defined as 



X'N = Ni 



Atv - (1 - y/CNT^ 



-(1 - Vc^)3C^ 

The above test is particularly suited to the model of Section 
III-BI in which the matrices Pk are non-positive definite when 
fik = and Qffc = — 1 for all k, corresponding to a sudden 
drop of a zero mean random parameter d{k) to zero. 

When the matrices Pk have both positive and negative 
eigenvalues, a deterministic choice has to be made by the 
experimenter. In the most general setting, to ensure a false 

'We recall that the false alarm rate is the probability of declaring IKq under 
true hypothesis JCq. 



alarm rate lower than 77, one has to choose two scalars 

a{r]), hijj) G M U {—00, 00} such that 



p({A;>a(77)}u{A^>5(r;)}) < r;. 



(15) 



The choice of 0(77), 6(7/) depends primarily on the structure of 
Pk and will impact the correct detection rate for fixed false 
alarm rates. 

In ll40l . the asymptotic independence of the fluctuations of 
the largest and the smallest eigenvalues of GUE matrices is 
proved, while the same result for the eigenvalues Ai and Aat 
of EE* under IKq is conjectured. Following this conjecture, 
( fTSl ) would become asymptotically 

T2{a{ii))T2{b{r,))>l-ij. 

For any fixed h, taking b{rj) = b, the hypothesis test now 
becomes 

min U - A'^, (T2)-i (l—^) ~ X[} t 0. 



T2ib) 



JCo 



In particular, for b — 00, T2{b) = 1 and then the test reduces 
to 

which is the same test as proposed in ( fTsT l. Taking instead 
a{ri) = —00, we obtain the test (fT4l i. 

For rather symmetrical distributions of the eigenvalues u!k,i 
of Pk around zero, it may be interesting to set b{ri) — a{r]), 
in which case 



b{v) - (^2)-^ (v/T^) 



In this setting, the decision test is now 

J-Cq 






In the following section, we assume that the procedure of 
failure detection was achieved successfully and that we are 
now interested in localizing the failures. 



(14) B. Localization algorithm 

We now wish to detect all possible failure events from a 
set of failures indexed by fc e {1,...,K}. The index set 
{1,...,K} may gather all events accounting for a single, 
as well as multiple, local failures. Similar to the previ- 
ous sections, we denote pk^i == 1 + ujk.i + c — ^^^ and 



Ck.i = -, ^, we define the mapping "Kk to be such 

that OCkii) = jk,i + ••• + jk,i-i if 1 < J < Sfc and 
Xkii) = N- [jk., + ...+ jk.tu) if Sfe + 1 < ^ < tfe. Finally, 
we denote life ^ any projector on the subspace generated by 



3Cfc(i)+l,---,AgCfe(i)+ji, 



the eigenvalues A 

An initial hypothesis rejection may be performed at this 
stage to select only those hypotheses "Kk such that the pk.i 
are consistent with the observations Ai, . . . , \n- For instance, 
if the largest eigenvalue Ai is significant in the system model. 



11 



one may preselect from the set {1, , 
indexes defined by 



arg ^ 



(fei 



mm 

,fcl,)C{l,, 



max 

,K} ki,...,ki 



, K} the L hypothesis 



1^1 - Pkj,l\ 



However, if K is large, many hypotheses may have very 
close parameters p^ ,, so that it is hazardous to conclude on 
the most likely hypothesis based only on the eigenvalues of 
EE*. However, since different Pj. matrices have in general 
very distinct eigenspaces, we propose the following subspace 
localization test, which decides on the hypothesis IK^* for 
which fc* is given by 



k* = argmaxgfc (^{V^n,^^^)^ 



(16) 






with gk the actual density of the vector 

^{Pk-, gfc) = {1, . . . ,Pfc, gfc, • ■ • , Tk}, S the set of (remaining) 
indexes k such that L{pk,qk) is non-empty, and where 

rk 



y 7 



NUt 



k.i 



(n. 



Ck,iljk.i Uk 



and 



7-fc A 



AgCfc(i)+i - Pk,i 



Xn 



''^k{i)+Jk.i - Pk.> 

Note that we need here to specify the indexation i E L{pk,qk) 
since we do not assume A5. 

From Theorem |4] this probability can be approximated for 
large n, which provides immediately a maximum likelihood 
test for the most asymptotically likely 'Kk hypothesis. In 
the particular case where the Uk^t all have multiplicity one, 
according to Corollary ^ as N,n grow large, the vectors in 
the test (fTSI l are asymptotically independent and Gaussian. We 
therefore substitute the test (fTSI l by the following test, leading 
to the estimator k defined as 

n f{iV^':n,LlJ■,C{pk,^)) (17) 

ie£(pfc,gfc) 

where f{x;V,), x G C™, il G C™^™, is the m-vai-iate real 
normal density of zero mean and covariance C at point x, and 
C{pk,i) is defined similarly as in Corollary |4] with uji replaced 
by uJk,i- Taking a log-likelihood notation, this is explicitly 



k = arg max 

kes 



k = arg mm > 

kes ^-^ 

l£C{pk,qk) 



+ logdctC(pfe,,) + 21og(27r)]. 



In the general case where ujk,i has multiplicity jk^i, as 
discussed previously, it is simpler to restrict the detection 
test to the eigenspace projections (^''„)i(E£(pfc.<;fc) only. In 

this case, denoting yA„ = {V^^^n,D^'^y^^,n,u) ^ l^^"'" with 
G R-""' the vector of the diagonal entries of V^^ and 
jk,i) jjjg vector of the real and imaginary parts 
of the upper-diagonal entries of V^^^, we obtain the test 

k = arg mill 



T/fe 



fees ^-^ I 
ieL(pk,qk} 



V ' i.nJ * i.n 



[C{pk,^)]ll 

HLlogi[C{pk,.) 



' Jk 



2 log(2^) 



(18) 



Remark 3: We provide below some remarks and discuss the 
advantages of the detection tests proposed in Section lTV-AI and 
the localization algorithms (fTTIi-lfTSll compared to the optimum 
maximum likelihood approach: 

• the detection algorithms proposed in Section IIV-AI are 
very versatile, as they adapt to multiple failure scenar- 
ios showing small rank perturbations in the population 
covariance, and provide a theoretical expression of the 
minimum ratio cn — N/n necessary for detectability; 

• unlike the traditional maximum-likelihood approach 
which tests the joint distribution of S for all hypotheses 
\ < k < K, and therefore leads to calculus of the order 
0{N'^) (or O(iV^) with some simplification methods) for 
each fc, the proposed localization algorithm ( fTTb is based 
on a test requiring for each k (taken from a possibly 
reduced subset of {1, ... , K}) eigenvector projections of 
computational load of order 0{N); 

• we may decide not to consider the joint fluctuations 
of all eigenvalues found outside S-^, but only some 
of them. This leads to an asymptotically less efficient, 
although much faster, algorithm, where Z {pk ,qk) in ( flTb 
is replaced by L{p',q') for given p' < pk, q' < qk, for 
all k. For N not too large, it is in fact preferable to con- 
sider only a few eigenvalues and eigenspace projections 
simultaneously, due to convergence speed limitations of 
the limiting normal distributions; 

• the entries L*^„ of the vector (yj*j^,L*^„) may also be 
discarded, especially in scenarios where eigenvalues of 
Pk are very similar for each hypothesis Jik- This may 
again increase the convergence speed of the asymptotic 
approximation for not-too-large N, while it is expected 
to perform worse for large N. 

So far, we have performed failure detection under the 
important assumption that the failure scenarios form a discrete 
set {1, . . . ,K}. This assumes in particular that the failure 
amplitudes are known prior to detection and localization. In 
the next section, we use Corollary |4] to improve this approach 
in the particularly simple example of Section III-BI when the 
failure amplitude is a priori unknown. 

C. Extension to unknown failure amplitude 

In this section, we assume the scenario where the eigen- 
vectors of the perturbation matrix Pk are independent of the 
amplitude of the failure parameters, in the sense that a change 
in magnitude of the failure of type k does not affect the 
eigenspaces of Pk- This is for instance the case of the single- 
failure scenario of Section III-BI for which we recall that 
Pk expresses as Pk = /3kR'~^ HekelH* R^^ with /3k the 
failure parameter We now assume /3k unknown, which is a 
more realistic assumption than assuming it perfectly known in 
advance. We also suppose that X has i.i.d. Gaussian entries. 
Based on a simple extension of the algorithm presented in 
Section IIV-BI to unknown ojk, we provide hereafter a second 
localization algorithm. 

For notational convenience, we assume Pk — ujkUkul for 
each k and that ujk > \fc, unknown. We then denote A the 
largest eigenvalue of EE* and u its associated eigenvector. 



12 



(3.31) (4.25) (4.29) (4.41) (2.82) 

'IV 0.80 -C^y- 0.77 -(eV 0.84 -(J)- 0.96 -(^ 



A A 

0.78 0.92 0.81 0.78 0.85 0.96 0.82 0.99 0.99 

(ij- 0.91 -Q- 0.83 -(?)- 0.95 -{t)- 0.96 -(^ 
(2.36) (4.50) (4.12) (4.43) (3.71) 

Fig. 2. Network of Af = 10 sensors. The correlation E[y(i)*j/(j)] between 
data on sensors i and j, i 7^ j, can be read on tlie link (j, j), while E[|?/(i)p] 
variances are shown in parentheses. 



Obviously, since ujk is not known, neither is pk ■ Therefore, 
we cannot proceed here to localization based on the fluctua- 
tions of A. Instead, we will use A precisely as an estimate of 
Pk, which we know is consistent with growing N, n. From A, 
assumed larger than (1 + y/c)'^, we want to derive an estimate 
uj of ujk {k is the effective failure index). This is obtained from 
an inversion of the relation (|7|i. Precisely, we obtain 



co^l{\-{l + c)) + y{\-il + c)r-4c 



if A > (1 + Vc)^ and 



if A < (1 - V^)2. 

From this estimate, we then obtain an estimate ( of (k as 
follows 

1 - CW-2 



C 



1 + CUJ 



A natural object to consider for the failure localization is 
now JM^up — C- To provide a diagnosis test, we need to derive 
the fluctuations of this random variable. From Theorem H] 
the fluctuations of \/]V(|u^.up — Cfc) depend on Uk but not 
on Ufc. From the expression of (^, it is immediate that the 
fluctuations of vN{C — (k) also depend on ujk only. But 
since ujk is estimated by a), irrespective of the failure index k, 
the diagnosis test leads to finding the most likely argument 
k' among K variables with same Gaussian statistics. This 
therefore simplifies the estimator k' of the most likely index 
k to the following minimum-distance estimator 



k' — arg inin 

ke{l,...,K} 



-c 



Note importantly that, contrary to our proposed scheme, 
the optimal maximum-likelihood localization method cannot 
be easily extended to the scenario of unknown failure ampli- 
tude, therefore bringing another significant advantage of the 
subspace approach. 

In the next section, we provide simulation results for 
single failure localization for the detection and localization 
algorithms assuming the failure amplitude known or unknown, 
applied to the scenarios of Section III-AI and Section III-BI 
respectively. 



D. Simulations 

In this section, we focus on the application of the algorithms 
designed in Sections IIV-AI and IIV-BI for single node failure in 
the scenario of Section III-AI and single parameter change in 
the scenario of Section III-BI 

1) Node failure in the scenario of Section 1/7- A 1 Our first 
application example relates to the sensor network model 
y ^ HO + aw of Section III-AI for A^ = 10 nodes, p = N, 
and a'^ = —20 dB. This is depicted in Figure |2] where 
the entries of HH* + ct^/at are presented. We also take 
aj. = ^i^i{HH*)ki, which is a natural assumption to avoid 
that a mere energy detector on y{k) provide a simpler solution 
to our problem. This failure amplitude is assumed known 
by the experimenter. In practical scenarios, this may arise 
if a sensor starts returning time delayed data, supposedly 
uncorrected with real-time data but with same variance. We 
assume a single failure scenario. In this context, it appears 
that, for all k, cuk^i > 0, uJk.2 < and ujk.i is much 
larger than |a;fc,2|- It is therefore more interesting only to 
consider the largest eigenvalue of SS* to detect and locate 
an hypothetical node failure. Under these conditions, the 
theoretical threshold for cn = N/n (if N, n were large) is 0.8 
with the worst-case failure corresponding to a failure of node 
10. We therefore carry out 100 000 Monte Carlo simulations 
of node 10 failures for n varying from 8 to 140 and under 
false alarm rates varying from 10^^ to 10^"*. This is depicted 
in Figure [51 where it can be observed that, for n = 8, detection 
and localization are barely possible, although it is clearly 
the starting point where detection becomes feasible. For not 
too large n, while detection rates increase, we observe that 
localization capabilities are still unsatisfying. This is mainly 
due to the inappropriate fit of the large dimensional model 
with A^ = 10 and with the eigenvectors corresponding to 
the extreme eigenvalues of SS* being too loosely correlated 
to their associated population eigenvectors. Larger values 
of n show much better performance with miss localization 
probability going to zero as n —>■ oo. In particular, about 
five times the asymptotically optimal ratio n/N is required 
for localization to be very efficient. In this case, the large 
dimensional model for the fluctuations of the eigenvalues and 
eigenvectors is more adapted. 

The same conditions are simulated for a system with N = 
100 nodes in which each node has eight neighbors and with 
correlation values of the same order of magnitude as in Figure 
|2] The detectability threshold for N/n is here 0.85 and we 
still consider the worst case failure scenario. This is depicted 
in Figure HI where one can see that smaller ratios n/N over 
the asymptotically optimal threshold are demanded for high 
detectability and localization ability to appear, when compared 
to the scenario A = 10. 

2) Sudden unknown parameter change: In this section, we 
consider the parameter change scenario of Section III-BI We 
still consider the network of Figure |2] and a^ = —20 dB 
as above. We now assume a sudden change of parameter 
0(10) with /3io = 2, being the worst case scenario for failure 
identification if f3k = 2 for all k. We depict the performance of 
the failure detection and locaUzation algorithms and compare 
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Fig. 3. Correct detection (CDR) and localization (CLR) rates for different 
levels of false alarm rates (FAR) and different values of n, for node 10 failure 
in the sensor network of Figure |2] The minimal theoretical n for observability 
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Fig. 4. Correct detection (CDR) and localization (CLR) rates for different 
levels of false alarm rates (FAR) and different values of n, for worst case 
node failure in a 100-node sensor network. The minimal theoretical n for 
observability is n = 85. 



the settings where fik is known or unknown in advance to the 
experimenter. In the former scenario, we apply the locaHzation 
algorithm of Section IIV-BI based on the joint fluctuations of 
the extreme eigenvalues and eigenspace projections, while in 
the latter, we apply the localization algorithm of Section ITV-CI 
where a prior step of eigenvalue inference is performed before 
the study of the fluctuations of the eigenspace projections. The 
results are presented in Figure |5] 

It appears from Figure |5] that the suboptimal algorithm of 
Section HV-CI performs only slightly worse than the algorithm 
of Section lTV-Bl for large n, and that it even performs better for 
small n. This last observation is explained by the inadequacy 
of the theoretical value of Cfe for too small values of n. It 
is therefore interesting to see that, for practical purposes, the 
absence of prior knowledge on the amplitude of the failure 



Fig. 5. Correct detection (CDR) and localization rates for different levels 
of false alarm rates (FAR) and different values of n, for sudden change 
of parameter 10 in the Scenario of Section III-BI Comparison is made 
between localization assuming /3fc known (CLR) and localization assuming 
/3fe unknown (CLR-2). The minimal theoretical n for observability is n = 22. 



does not severely reduce the efficiency of the localization 
algorithm. 



V. Conclusion 

In this article, a characterization of the joint fluctuations of 
the extreme eigenvalues and corresponding eigenspace projec- 
tions of a certain class of random matrices is provided. This 
characterization was used to perform fast and computationally 
reasonable detection and localization of multiple failures in 
large sensor networks through a general hypothesis testing 
framework. The main practical outcomes of this article lie first 
in a characterization of the minimum number of observations 
necessary to ensure failure detectability in large networks and 
second in the design of flexible but simple algorithms that can 
be adapted to multiple types of failure scenarios consistent 
with the small rank perturbation random matrix model. We 
also extend the detection and diagnosis approach to scenarios 
where the amplitudes of the hypothetical failures are not a 
priori known. Practical simulations suggest that the proposed 
algorithms allow for high failure detection and localization 
performance even for networks of small sizes, although for 
those much more observations than theoretically predicted are 
in general demanded. 



Appendix A 
Proofs of results of SectionHIII 

A. Proof of Lemma |2] 

We shall assume without loss of generality that i = 1. In 
Section IIII-B2I we saw that 

f/rfiiC/i = ^ / MzyH{z)-^Mz)dz 



14 



with probability one, where 

li(z)* =zU*ilN+Py^Qiz)U 
-UlQ{z)U 



(1+Wl)2 



and 



A2{z) = n{Ir + ^)-^U*Q{z){In + Py-^Ui 

= (1 + wi)-^f](/^ + n)-'^u*Q{z)Ui 

(take 61 and 62 as any two columns of Ui in ([8])). Similarly, 
1 



where 






Ai{z)* ^z'm{z)Ul{I + P)-^U 
-T^Ah, 0], 

(1 +L0l)2 

A2{z) = m{z)Q{i + ny^u* {I + py^Ui 



u}im{z) 



(l+Wi)3/2 

Fixing z G 71, we have 

Hiz)-' = H{z)-^ - H{z)-^E{z)H{z)-^ + 0{\\E{z)f) 
where 

E{z) = ^(z) - H{z) = zVt{Ir + VL)-^U*{Q{z) - m{z)I)U. 
Guided by this observation, we write 



Vi,n 



N 
'2m 



71 



(Il(z)*ij(z)-1A2(Z) 



Ai{z)*H{z)-^A2{z)\ dz (a.s., n large) 



= "^/ iM^y - M^r)H{z)-^A2{z)dz 

+ ^/ Ai(z)*iJ(z)-i(l2W-^2W)d;^ 

-^ / Ai(z)*iJ(z)-i£;(z)if(z)-iA2(z)dz 
27ri /^^ 

+ £ 

= Zi + Z2 + Z3 + £ 

where £ contains all the higher order terms that appear when 
we develop the integrand at the right hand side of the first 
equality. 

In what follows, we successively study each of the terms 
at the right hand side of this equation. Recalling ^, the term 
Zi writes 



the residue theorem and the identity (1 + loi) 
h{pi))/h{pi), we obtain 

_ i + Mpi) 

Similarly, the term Z2 writes 



NUl{Q{pi) - m{pi)I)Ui 



Z2 — Zi — 



l + h{p^) 



h'ipi) 
Turning to Z3, we have 



iVf7*(Q(pi)-m(pi)/)[/i. 



z2m(z)2 U*{Q{z)-m{z))Ui 



rdz 



27ri J^^ (1 + uji) ((1 + uji)/uJi + zm{z)y 

which shows that we have a pole with degree 2 in Int(7i). 
Write the integrand as G{z)/g{z) and recall that the residue of 
a meromorphic function f{z) associated with a degree 2 pole 

at zq is \imz^zad[{z — zo)'^f{z)) /dz. After some simple 
calculations, this results in 



Z.= 



G{p,)g"{p,) G'{p,) 

9'{pif 9'{piY 

l + /i(pi) fh{pi)h"{pi) 



-2 



h'ip,) V h'ip.r 
X ^/NUliQipi) ~ m{pi)I)Ui 
/i(pi)(l + /i(pi)) 



^,(^^^2 .NUt{Q'{p,)-m'{p,)I)U,. 

We now show briefly that the last term £ in the expression 
of Vi „ converges to zero in probability. A more detailed 
argument is given in |29|. Recall that £ accounts for all the 
higher order terms that show up when we expand the integrand 
AiH^^A2 — AiH^^A2. Let us focus on one of these terms, 
namely 



H{z)-^ - H{z)-^ + H{z)-^E{z)H{z)-^\ A2{z)dz 

and show that it converges in probability to zero. The other 
terms can be treated similarly. First, we can show that 

\\H{z)-^ - H{z)-^ + H{z)-^E{z)H{z)-^\\ < K\\E{z)f 

on 71 where K is some constant. Now we write 

E{z) = zn{Ir + ny^U*{Q{z) - a{z)I)U 
+ z{a{z) - m{z))n{Ir + ny^ 
= Ei{z) + E2{z). 

Noticing that Ai{z) and ^2(2;) are bounded on 71, and writing 
z — pi + Rexp{2nTd) on 71, the result is shown if we show 
that ^ 



\E,{pi+Re' 







(19) 







Zi 



N f zm,{z) Ul{Q{z) - m{z)I)Ui 
2^ J J, il+LOi)iil+LOi)/uji + zm{z)) 



dz. 



The denominator has one simple zero in Int(7i). With prob- 
ability one, the numerator has no zero in Int(7i). Using 



for i = 1,2. Lemma [U shows that E||i?i(z)|p < K'/n on 
71 where the constant K' is independent of z. By Markov's 
inequality, ( fT9l ) is true for Ei . Convergence for E2 is obtained 
from Assumption A4 in conjunction with the analyticity of 
a{z) — m{z), as shown in |29|. 

Taking the sum Zi + Z2 + Zs, we obtain the desired result. 
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. Proof of Lemma |5] 






Set i = 1 and write 


n{ir + n)-^ = 


^1. 


et y„ = pi + xj^/N . Write U = 


C/i C/i 



ii^(2/n) 



and 

In + f^llnU^QiVnWl /3iynU*Qiyn)Ui ^ 

VnBiUlQ{yn)Ui Ir-n+VnBiUtQ{Vn)Ui 

H21 H22 
Let us study 

N^ det i?(y„) = det(i722) det{VNHn - ^Hi2H^2^H2i). 
By Lemma [T] 

"(1 - ^2/Pi)In 



Ho 



hence 



(1 - A//3i)/,,j 
AciH22^X{{l-pilPiy'. 



i>i 



By the same lemma. 



E 



Nl,^XX')e[a-e^b+e]\\Hl2\\^ <K/VN 



and 



E 



UxX')e[a-e^b+e]\\H2l\\\ <K'/VN 

where K and K' are some constants, hence 

det{VNHii - VnHi2H22^H2i) - det(V^^ii 
Let us study this last term. We write 

V^vifii - V7V/3iy„C/r(Q(yn) - a{yn)lN)Ui 
+ y/N I3i{yna{yn) - pia{pi))Ij, 
+ VNl3i{pia{pi) - pim{pi))Ij, 
= ^1 + Z2 + ^3. 
Writing 

Qivn) - g(pi) = Q(y„)(Q(pi)-^ - Qiynr')Qipi) 

we have 

Zi - %/lV/3iy„[/*(g(pi) - aipi)I)Ui 
= a:/?iy„f/*(g(2/„)Q(pi) - 7V-i(trQ(y„)Q(Pi))/)C/i 

by Lemma [T] From A4, we further have 

Zi - %/iV/3ipif/*(Q(pi) - m(pi)/)f/i A 0. 
Turning to the second term, we can show using A4 that 
yna{yn) - pia{pi) 



0. 



Pix- 



x/VN 



f3ix{pim{pi)y . 



Again by A4, Z3 ^-4 0. This results in 



e>i 



X det [VNPipiU*{Q{pi) - m{pi)I)Ui + Pixh'ip,) 
The same argument for i > 1 leads to the result. 



C. Proof of Lemma |4] 

Recall that XX* admits the spectral factorization XX* = 
WAW* where W and A = diag(Ai, . . . , Aat) are indepen- 
dent, and where W is Haar distributed on the group of N x N 
unitary matrices. 

From Assumption A4 and the analyticity of /i, we can show 
as in |,36J that 



\ k=l •' 



/i(A)d^(A) 







hence the lemma is shown if we show the result on 

'u*f,{XX*)U, - U Eti /.(Afe)l /,.' 
U*g,iXX*)U, - h Ef=i g.iXk)] In 



S.-n 



Let Z ht a. N X r random matrix with independent CINf(0, 1) 
elements, chosen independently of A. Write Z ~ [Zi . . . Zt] 
where the blocks Zg have the same dimensions as the Ue. Then 
S_^ is equal in distribution to 



N 



[Z{Z*Z)- 
[Z{Z*Z)- 



V.(A) - 'lM^\ {Z{Z*Z)-% 



where \Z{Z* Z)^^\i is the matrix formed by the columns ji + 
. . . + ji-i to ji + . . . + ji of Z{Z*ZY^. By the law of large 
numbers, N^^Z*Z ^-4 7^, hence it will be enough to show 
the result on 

o _^, Z*{f,{A)-N-HYfdA)lN)zy^' 
^' Z*(g^,{A)-N-Hrg,iA)lN)z, 



N 

E 
fe=i 




'^i,k^i,k 



/ J / t — 1 

Zi.n]- Write Ci,fc = 



where we have written Z* = [zi^i 

/,(Afe) - iV-itr/,(A) and d^,k = 9^{Xk) - N-Hvg,{A) 

Up to an element rearrangement, Sn can be rewritten as the 



1 ^ 



Vih ' 



Ci,k 

di.k 



(20) 



where for every i — 1, . . . ,t, the N vectors Vi k are valued in 
Wi . We shall show that this sum converges in law to a Gaus- 
sian R^(-'i"' ''■'*) -valued random vector whose covariance 
matrix is equal to the covariance matrix of ( fT2b rearranged 
similarly to Sn- 

We observe that the summands of ( l20l) are centered and are 
independent conditionally to A. Observe also that for every k, 
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the vectors {vi^k)i=i are independent and that the elements of 
each of these vectors are decorrelated. Based on A2 and A3, 
we have 



1 ^ 



k=l 



diae; L, 



- 


Vi,k <E) 


Ct.k 




t / 


Wj,fe ® 


'Ci^k 

_di,k 




t \T 


A 



N r 



AtZ^ 



k=l 






i=l 



diag (/j^ ® i?i)*=i 



which coincides with the covariance matrix of (fT2l) after the 
rearrangement. 

Furthermore, thanks to A3, it is easy to see that the 
Lyapunov condition 



N 



1 ^ 



fc=i 



- 


Vi,k ® 


Ct,k 

dt,k 




t 


2+2J7 


A 



is valid for any 77 > 0, hence ( |20] | satisfies the conditions of 
the central Umit theorem, which proves the lemma. 
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